library(tidyverse)
library(knitr)
library(plotly) ; library(viridis) ; library(gridExtra) ; library(RColorBrewer) ; library(ggpubr)
library(biomaRt)
library(polycor)
library(caret) ; library(ROCR) ; library(car) ; library(MLmetrics)
library(corrplot)
library(expss) ; library(knitr) ; library(kableExtra)
library(foreach) ; library(doParallel)

SFARI_colour_hue = function(r) {
  pal = c('#FF7631','#FFB100','#E8E328','#8CC83F','#62CCA6','#59B9C9','#b3b3b3','#808080','gray','#d9d9d9')[r]
}


Load Dataset


clustering_selected = 'DynamicHybrid'
clusterings = read_csv('./../Data/clusters.csv')
clusterings$Module = clusterings[,clustering_selected] %>% data.frame %>% unlist %>% unname
assigned_module = clusterings %>% dplyr::select(ID, Module)

# Dataset created with 20_04_07_create_dataset.html
dataset = read.csv(paste0('./../Data/dataset_', clustering_selected, '.csv'), row.names=1)
rownames_dataset = rownames(dataset)
dataset = dataset %>% mutate(Module = clusterings$Module, gene.score = as.character(gene.score)) %>%
                      mutate(gene.score = ifelse(gene.score=='Others', 'None', gene.score)) %>%
          dplyr::select(-matches(clustering_selected))
rownames(dataset) = rownames_dataset

# Fix any Gene Significance that is NA
GS_missing = rownames(dataset)[is.na(dataset$GS)]
if(length(GS_missing)>0){
  # Gandal dataset
  load('./../Data/preprocessed_data.RData')
  datExpr = datExpr %>% data.frame
  
  for(g in GS_missing){
    dataset$GS[rownames(dataset) == g] = polyserial(as.numeric(datExpr[g,]), datMeta$Diagnosis)
  }
}



# SFARI Genes
SFARI_genes = read_csv('./../../../SFARI/Data/SFARI_genes_01-03-2020_w_ensembl_IDs.csv')
SFARI_genes = SFARI_genes[!duplicated(SFARI_genes$ID) & !is.na(SFARI_genes$ID),]


# GO Neuronal annotations: regex 'neuron' in GO functional annotations and label the genes that make a match as neuronal
GO_annotations = read.csv('./../Data/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>% 
              mutate('ID'=as.character(ensembl_gene_id)) %>% 
              dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
              mutate('Neuronal'=1)


# Add gene symbol
getinfo = c('ensembl_gene_id','external_gene_id')
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl',
               host='feb2014.archive.ensembl.org') ## Gencode v19
gene_names = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), values=rownames(dataset), mart=mart)


rm(getinfo, mart, rownames_dataset, GO_annotations, g, GS_missing)


The features that will be considered for the classification model will be the ones WGCNA uses to identify significant modules and genes:



Filtering the 5530 genes that were not assigned to any cluster (represented as the gray cluster)

rm_cluster = dataset[is.na(dataset$MTcor),'Module'] %>% unique %>% as.character

new_dataset = dataset %>% filter(Module != 'gray' & !is.na(MTcor)) %>% 
              dplyr::select(-c(matches(paste('pval|Module')), MMgray)) %>%
              mutate('absGS'=abs(GS), 'SFARI'=gene.score!='None') %>% dplyr::select(-gene.score)
rownames(new_dataset) = rownames(dataset)[dataset$Module != 'gray']

rm(rm_cluster)


Summary of the changes made to the original WGCNA variables:


  • Using Module Membership variables instead of binary module membership

  • Including a new variable with the absolute value of GS

  • Removing genes assigned to the gray module (unclassified genes)

  • Adding the Objective variable: Binary label indicating if it’s in the SFARI dataset or not

table_info = new_dataset %>% head(5) %>% t 

table_info %>% kable(caption = '(Transposed) features and their values for the first rows of dataset', 
                     col.names = colnames(table_info)) %>% kable_styling(full_width = F)
(Transposed) features and their values for the first rows of dataset
ENSG00000005243 ENSG00000143815 ENSG00000171100 ENSG00000104783 ENSG00000163288
MTcor 0.5467749 0.3344469 -0.0355657 0.3969497 0.3460181
GS 0.6059364 0.2858943 0.1872330 0.4001396 0.3973047
MM.00C19C -0.1206088 0.3654242 0.3787122 -0.0413128 -0.0507447
MM.00BBDA -0.2359082 0.1844031 0.1349984 -0.3245480 -0.0656331
MM.FF64B4 -0.4554077 -0.0006262 -0.0177737 -0.4268871 -0.1215246
MM.CD79FF 0.0810707 -0.0213142 0.0234191 -0.3155081 -0.0768211
MM.E58700 0.1362213 0.2644513 0.0013391 -0.2407284 0.1501250
MM.00BECE -0.1900722 -0.1956851 -0.1019865 -0.2720456 -0.0129917
MM.00BF7A -0.3865027 -0.3097687 -0.0594399 -0.4141432 -0.1767419
MM.F8766D -0.1296755 -0.2635687 -0.1777153 -0.2874802 -0.1709495
MM.00BFC8 -0.0553910 -0.1721149 -0.2508106 -0.5272773 -0.1668932
MM.F57962 -0.2945354 -0.2428728 -0.1361109 -0.6299000 -0.2289118
MM.00BD5B -0.0478729 -0.0820154 -0.1914338 -0.3941551 -0.0615747
MM.53B400 -0.1229521 -0.1243407 -0.2499512 -0.4366579 0.0167560
MM.BB81FF -0.1895681 -0.0444843 0.1995541 0.0350543 -0.1435966
MM.00BE70 -0.1898026 -0.1615888 -0.0400521 -0.2031248 -0.1765158
MM.689BFF -0.1349865 -0.0321037 0.1569906 -0.2934116 -0.2555161
MM.D49200 -0.0577336 0.0667018 0.1118843 0.2990581 0.1547791
MM.40B500 -0.0296442 0.1490564 0.1777671 -0.3071037 0.1905981
MM.9DA700 -0.0264936 -0.0590357 0.1583884 -0.2352978 0.1241876
MM.00B4F0 -0.1321228 -0.2232373 -0.3451113 -0.1545722 -0.0133171
MM.ABA300 -0.3497061 -0.2139791 -0.2626484 -0.3176980 -0.1013731
MM.24B700 -0.1780637 0.1614555 0.2574047 -0.1497984 -0.1423753
MM.EC823B 0.0414052 -0.0415483 0.0416945 0.1367569 -0.0430319
MM.00BDD4 -0.0772010 0.0964134 0.1092253 -0.1698598 -0.0251578
MM.8CAB00 -0.2652459 -0.0374433 -0.0481687 -0.3364900 -0.1422723
MM.00B81C -0.1159573 0.0459810 0.0941499 0.1305122 0.0134355
MM.00C083 -0.2088419 -0.0753278 -0.0235310 -0.1306668 0.0312405
MM.95A900 -0.1835726 -0.0217957 -0.1836436 -0.1239786 -0.0261126
MM.A4A500 -0.2212386 -0.1894831 -0.1860289 -0.3889673 0.1099064
MM.FF62BB -0.1403988 -0.0934565 -0.0927384 -0.3104293 0.0019839
MM.00C08B -0.0192068 -0.1675827 -0.0985837 -0.4207551 -0.0889154
MM.00C0BA -0.2078071 -0.3823642 -0.1488761 -0.4712894 -0.2585668
MM.00C0B3 -0.4456070 -0.1253304 0.0406955 -0.5450825 -0.3430102
MM.B2A100 -0.3446382 -0.1240596 -0.0454376 -0.4892012 -0.1666268
MM.7B97FF -0.2555002 -0.1967161 -0.2585973 -0.2979366 -0.0174912
MM.00BA42 -0.2729022 -0.2981534 -0.0916584 -0.3625126 -0.1905154
MM.E18A00 -0.4274569 -0.4073737 -0.2346541 -0.3578283 -0.2519161
MM.00B8E6 -0.3366292 -0.1700978 0.1198042 -0.0816026 -0.1078875
MM.61B200 -0.4020531 -0.1155759 -0.1045518 -0.1189860 -0.0682339
MM.D575FE -0.2349892 -0.2093592 -0.2886493 -0.0936652 0.0445961
MM.00BFC1 -0.2367516 -0.3145291 -0.0806798 0.0715188 -0.0895351
MM.2FA2FF -0.3798354 -0.3264221 -0.2342503 -0.0671347 -0.1496347
MM.E36EF6 -0.1875503 -0.3140448 -0.1320553 -0.0669857 0.1356024
MM.FF699B -0.6363982 -0.1026271 0.0214001 -0.4129611 -0.2573844
MM.519EFF -0.4208855 -0.3441752 -0.2570868 -0.2611152 -0.0156430
MM.998EFF -0.3201243 -0.2779016 -0.3125953 -0.4686434 -0.0520988
MM.E98428 0.1531579 -0.0217965 0.0207042 -0.0745673 -0.0930563
MM.FA7377 0.1197057 -0.0277990 0.0893684 -0.0120158 -0.2487928
MM.FE6E8A -0.0906448 -0.0184753 -0.0258878 -0.2399785 -0.3425219
MM.79AF00 0.0868363 -0.1747982 0.0176828 -0.2155793 -0.1177606
MM.FF65AC -0.1189480 -0.1360682 -0.0694710 -0.3204034 -0.1300504
MM.00A6FF 0.0403894 -0.1802565 -0.2234705 -0.0282003 -0.2103655
MM.00A9FF -0.2015416 -0.3060135 -0.2758527 -0.3706763 -0.2408232
MM.D88F00 -0.2100707 -0.1914377 -0.1657217 -0.3101389 -0.1473895
MM.FC7181 -0.0358048 -0.2969291 -0.2665994 -0.1912806 -0.1325468
MM.A58AFF -0.0909786 -0.2148038 -0.2216710 0.0051962 -0.1903727
MM.F37C57 0.0995576 -0.1305354 -0.1050885 0.0493699 -0.0697874
MM.00B6EB 0.1548395 -0.1193182 -0.2650779 -0.0380516 -0.0142380
MM.83AD00 0.1676396 -0.0349410 -0.0100732 -0.1481515 -0.1198566
MM.00AFF8 -0.0754264 -0.0201151 -0.1152253 -0.1658506 -0.3036295
MM.00C094 0.1150228 -0.1177833 0.0558893 -0.0390213 -0.1793569
MM.DD8D00 -0.0669987 -0.1403359 -0.1984024 -0.1777784 -0.2149976
MM.F962DD 0.2843169 0.2954299 0.0437415 0.4454934 0.0893491
MM.00ACFC 0.1885257 0.1088179 0.0003510 0.1129265 -0.0533085
MM.CF9500 0.1992014 0.1897098 0.0497466 0.1914940 -0.0007938
MM.FF61C9 0.2913960 0.0638250 -0.2564893 0.1077031 0.1655579
MM.DC71FA 0.3949472 0.2797443 0.0904335 0.2005073 -0.0077633
MM.C99700 0.1416387 -0.1390416 -0.2812743 0.0247617 -0.0807855
MM.F266E8 0.3568656 0.1380833 -0.0167623 0.0152389 0.0381506
MM.00BC50 0.1336532 -0.0461056 -0.2326611 -0.0840295 0.0563062
MM.B89E00 0.0967391 -0.0292006 0.2210237 -0.0454091 -0.0655840
MM.F564E3 0.2101397 0.0060840 0.2184868 -0.1196188 -0.1008739
MM.00BAE0 0.0563814 0.3295106 0.2361469 0.1588367 -0.0195894
MM.C57DFF 0.2183101 0.5451805 0.3479702 0.1460148 -0.0404573
MM.8B93FF 0.0202065 0.2426882 -0.0096652 -0.2069296 -0.0066285
MM.FB61D7 -0.3282170 0.1387925 0.2457945 -0.2430987 -0.3488064
MM.EF7F4A -0.0552639 -0.0548921 -0.0321955 0.1942000 0.0335942
MM.00BD66 0.1753791 0.0894774 0.0604607 -0.0675281 0.5070396
MM.E86BF2 0.3793787 0.2120334 0.1066273 0.1775474 0.5575256
MM.FD61D0 0.3301527 -0.0963991 -0.1499150 0.0901139 0.4669795
MM.00B933 0.1263057 -0.1065320 -0.2683421 0.0333771 0.4181965
MM.6EB000 0.4753561 0.1297254 -0.0563796 0.3045269 0.4875669
MM.00C1A4 0.5442029 0.2365080 0.1238710 0.5149762 0.2645018
MM.B186FF 0.4115817 0.1654591 0.1367287 0.5192951 0.1558556
MM.00C1AB 0.5511627 0.1439284 0.0224213 0.2575399 0.1658112
MM.C49A00 0.2868600 -0.0077317 0.0488221 0.0376683 -0.0065447
MM.BE9C00 0.4275323 0.0645554 0.0818880 0.5729344 0.2975774
MM.FF6C93 0.2365965 -0.1432882 0.1754847 0.2696715 0.1035714
MM.ED68ED 0.3679420 0.2722669 0.0267371 0.7939112 0.3314038
MM.FF61C2 0.1958136 0.0214576 -0.0824511 0.5318761 0.2083225
MM.00B2F4 0.5724066 0.1387165 -0.0715436 0.5209587 0.3605438
MM.FF67A4 0.5429656 0.2516709 0.0241876 0.6405827 0.3296417
absGS 0.6059364 0.2858943 0.1872330 0.4001396 0.3973047
SFARI 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
rm(table_info)
original_dataset = dataset
dataset = new_dataset

rm(new_dataset)

The final dataset contains 7632 observations (genes) and 95 variables



Exploratory Analysis


PCA of Variables


The Module Membership variables are grouped by Module-Trait correlation, with positive correlations on one side, negative on the other, and both SFARI and absGS are in the middle of both groups (just like with Gandal’s dataset)

pca = dataset %>% mutate(SFARI = as.numeric(SFARI)) %>% t %>% prcomp

plot_data = data.frame('ID'=colnames(dataset), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2],
                       type=ifelse(grepl('MM', colnames(dataset)),'ModMembership',
                            ifelse(grepl('SFARI', colnames(dataset)), 'SFARI',
                            ifelse(grepl('GS', colnames(dataset)), 'GS', 'MTcor'))))


mtcor_by_module = original_dataset %>% dplyr::select(Module, MTcor) %>% unique
colnames(mtcor_by_module) = c('ID','MTcor')

plot_data = mtcor_by_module %>% mutate(ID = gsub('#','MM.',ID)) %>% right_join(plot_data, by='ID')

ggplotly(plot_data %>% ggplot(aes(PC1, PC2, color=MTcor)) + geom_point(aes(id=ID)) +
         xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) +
         ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)')) +
         scale_colour_distiller(palette = 'RdBu', na.value = 'darkgrey') + theme_minimal() +
         ggtitle('PCA of variables coloured by Module-Diagnosis correlation'))
rm(mtcor_by_module, pca)


PCA of Samples


  • The two main patterns that seem to characterise the genes are their Gene Significance and the Module-Diagnosis correlation of their corresponding module

  • Mean Expression doesn’t seem to play an important role

  • SFARI Genes seem to be evenly distributed everywhere (perhaps they have a slightly higher distribution in the 2nd principal component?)

  • It’s not clear what the 2nd principal component is capturing

# Mean Expression data
load('./../Data/preprocessed_data.RData')
datExpr = datExpr %>% data.frame
mean_expr = data.frame('ID'=rownames(datExpr), 'meanExpr' = rowMeans(datExpr))

# PCA
pca = dataset %>% t %>% prcomp

plot_data = data.frame('ID'=rownames(dataset), 'PC1'=pca$rotation[,1], 'PC2'=pca$rotation[,2], 
                       'SFARI'=dataset$SFARI, 'MTcor'=dataset$MTcor, 'GS'=dataset$GS) %>%
            mutate(alpha=ifelse(SFARI, 0.7, 0.2)) %>% left_join(mean_expr, by='ID')

p1 = plot_data %>% ggplot(aes(PC1, PC2, color=MTcor)) + geom_point(alpha=0.4) + scale_color_viridis() + 
     theme_minimal() + ggtitle('Genes coloured by Module-Diagnosis correlation') +
     xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1]),'%)')) +
     ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2]),'%)')) +
     theme(legend.position='bottom')

p2 = plot_data %>% ggplot(aes(PC1, PC2, color=GS)) + geom_point(alpha=0.4) + scale_color_viridis() + 
     xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1]),'%)')) +
     ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2]),'%)')) +
     theme_minimal() + ggtitle('Genes coloured by Gene Significance') + theme(legend.position='bottom')

p3 = plot_data %>% ggplot(aes(PC1, PC2, color=SFARI)) + geom_point(alpha = plot_data$alpha) +
     xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1]),'%)')) +
     ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2]),'%)')) +
     theme_minimal() + ggtitle('Genes coloured by SFARI label') + theme(legend.position='bottom')
p3 = ggExtra::ggMarginal(p3, type='density', groupColour=TRUE, size=10)

p4 = plot_data %>% ggplot(aes(PC1, PC2, color=meanExpr)) + geom_point(alpha=0.4) + scale_color_viridis() + 
     xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1]),'%)')) +
     ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2]),'%)')) +
     theme_minimal() + ggtitle('Genes coloured by mean level of expression') + theme(legend.position='bottom')

grid.arrange(p1, p2, p3, p4, nrow=2)

rm(pca, datExpr, datGenes, datMeta, dds, DE_info, mean_expr, p1, p2, p3, p4)



Dividing samples into Training and Test Sets


5.67% of the observations are positive. This can be a problem when training the classification model, so the samples in the training set should be balanced between classes before the model is trained.

table_info = dataset %>% apply_labels(SFARI = 'SFARI')

cro(table_info$SFARI)
 #Total 
 SFARI 
   FALSE  7199
   TRUE  433
   #Total cases  7632
rm(table_info)

To divide our samples into training and test sets:

set.seed(123)

sample_scores = dataset %>% mutate(ID = rownames(.)) %>% dplyr::select(ID) %>% 
                left_join(original_dataset %>% mutate(ID = rownames(.)) %>% dplyr::select(ID, gene.score), 
                          by = 'ID') %>% 
                mutate(gene.score = ifelse(is.na(gene.score), 'None', gene.score))

train_idx = createDataPartition(sample_scores$gene.score, p = 0.7, list = FALSE)
train_set = dataset[train_idx,]
test_set = dataset[-train_idx,]


rm(sample_scores, train_idx)


Label distribution in training set


To fix this class imbalance, we are going to use SMOTE, an over-sampling technique that over-samples the minority class (SFARI Genes) by creating synthetic examples, in the training set

cro(train_set$SFARI)
 #Total 
 train_set$SFARI 
   FALSE  5040
   TRUE  305
   #Total cases  5345


Labels distribution in test set


This set is used just to evaluate how well the model performs, so the class imbalance is not a problem here

cro(test_set$SFARI)
 #Total 
 test_set$SFARI 
   FALSE  2159
   TRUE  128
   #Total cases  2287


Logistic Regression


Train model

# https://shiring.github.io/machine_learning/2017/04/02/unbalanced
# https://topepo.github.io/caret/using-your-own-model-in-train.html#Illustration5


train_set = train_set %>% mutate(SFARI = ifelse(SFARI==TRUE, 'SFARI', 'not_SFARI') %>% as.factor)

k_fold = 10
cv_repeats = 5
smote_over_sampling = trainControl(method = 'repeatedcv', number = k_fold, repeats = cv_repeats, 
                                   verboseIter = FALSE, classProbs = TRUE, savePredictions = 'final', 
                                   summaryFunction = twoClassSummary, sampling = 'smote')

# Using ROC as metric because it doesn't depend on the threshold
fit = caret::train(SFARI ~ ., data = train_set, method = 'glm', family = 'binomial', metric = 'ROC',
                   trControl = smote_over_sampling)

# There is some perfect multicollinearity that doesn't let us do the vif analysis, so I'll remove those variables (you cannot use alias in caret::train, so I had to train the model again directly with glm)
ld.vars = attributes(alias(glm(SFARI~., data = train_set, family = 'binomial'))$Complete)$dimnames[[1]]

# Remove the linearly dependent variables variables
formula.new = as.formula( paste('SFARI ~ .', paste(ld.vars, collapse='-'), sep='-') )

# Retrain model without these variables
fit = caret::train(formula.new, data = train_set, method = 'glm', family = 'binomial', metric = 'ROC',
                   trControl = smote_over_sampling)


rm(smote_over_sampling, ld.vars, formula.new, k_fold, cv_repeats)


Performance


The model has an AUC of 0.6877

But the features are strongly correlated, which inflates the standard error of the coefficients, making them no longer interpretable, so perhaps it would be better to use another model

# VIF
plot_data = data.frame('Feature' = car::vif(fit$finalModel) %>% sort %>% names,
                       'VIF' = car::vif(fit$finalModel) %>% sort %>% unname) %>%
            mutate(outlier = VIF>10)

plot_data %>% ggplot(aes(reorder(Feature, -VIF), VIF, fill = !outlier)) + geom_bar(stat='identity') + 
              scale_y_log10() + geom_hline(yintercept = 10, color = 'gray', linetype = 'dashed') + 
              xlab('Model Features') + ggtitle('Variance Inflation Number for each Feature') + theme_minimal() +
              theme(legend.position = 'none', axis.text.x = element_text(angle = 90, hjust = 1))

rm(plot_data)

Correlation plot

corrplot.mixed(cor(train_set[,-ncol(train_set)]), lower = 'number', lower.col = 'gray', number.cex = .6, 
               tl.pos = 'l', tl.col = '#666666')


Possible solutions to Multicollinearity:


  1. Remove all variables with a VIF>10: We would lose all but two of our variables, not ideal

  2. Do Principal Component Regression: We would lose the relation between the prediction and the original features, which could be interesting to study

  3. Don’t do anything: Multicollinearity affects the coefficients and p-values of the regression, but it doesn’t affect the predictions, precision of the predictions or the goodness-of-fit statistics ref, but as with the previous option, we cannot study the coefficients of the regression

  4. Use Ridge Regression: The penalty it gives to high coefficients reduces the variance introduced by the correlation, making the coefficients interpretable again




Ridge Regression


Notes:

### DEFINE FUNCTIONS

create_train_test_sets = function(p, seed){
  
  # Get SFARI Score of all the samples so our train and test sets are balanced for each score
  sample_scores = dataset %>% mutate(ID = rownames(.)) %>% dplyr::select(ID) %>%
                  left_join(original_dataset %>% mutate(ID = rownames(.)) %>% dplyr::select(ID, gene.score), 
                             by = 'ID') %>% 
                  mutate(gene.score = ifelse(is.na(gene.score), 'None', gene.score))

  set.seed(seed)
  train_idx = createDataPartition(sample_scores$gene.score, p = p, list = FALSE)
  
  train_set = dataset[train_idx,]
  test_set = dataset[-train_idx,]
  
  return(list('train_set' = train_set, 'test_set' = test_set))
}



run_model = function(p, seed){
  
  # Create train and test sets
  train_test_sets = create_train_test_sets(p, seed)
  train_set = train_test_sets[['train_set']]
  test_set = train_test_sets[['test_set']]
  
  # Train Model
  train_set = train_set %>% mutate(SFARI = ifelse(SFARI==TRUE, 'SFARI', 'not_SFARI') %>% as.factor)
  lambda_seq = 10^seq(1, -4, by = -.1)
  set.seed(seed)
  k_fold = 10
  cv_repeats = 5
  smote_over_sampling = trainControl(method = 'repeatedcv', number = k_fold, repeats = cv_repeats,
                                     verboseIter = FALSE, classProbs = TRUE, savePredictions = 'final', 
                                     summaryFunction = twoClassSummary, sampling = 'smote')
  fit = train(SFARI ~., data = train_set, method = 'glmnet', trControl = smote_over_sampling, metric = 'ROC',
              tuneGrid = expand.grid(alpha = 0, lambda = lambda_seq))
  
  # Predict labels in test set
  predictions = fit %>% predict(test_set, type = 'prob')
  preds = data.frame('ID' = rownames(test_set), 'prob' = predictions$SFARI) %>% mutate(pred = prob>0.5)

  # Measure performance of the model
  acc = mean(test_set$SFARI==preds$pred)
  prec = Precision(test_set$SFARI %>% as.numeric, preds$pred %>% as.numeric, positive = '1')
  rec = Recall(test_set$SFARI %>% as.numeric, preds$pred %>% as.numeric, positive = '1')
  F1 = F1_Score(test_set$SFARI %>% as.numeric, preds$pred %>% as.numeric, positive = '1')
  pred_ROCR = prediction(preds$prob, test_set$SFARI)
  AUC = performance(pred_ROCR, measure='auc')@y.values[[1]]
  
  # Extract coefficients from features
  coefs = coef(fit$finalModel, fit$bestTune$lambda) %>% as.vector
  
  return(list('acc' = acc, 'prec' = prec, 'rec' = rec, 'F1' = F1, 'AUC' = AUC, 'preds' = preds, 'coefs' =coefs))
}


### RUN MODEL

# Parameters
p = 0.75

n_iter = 25
seeds = 123:(123+n_iter-1)

# Store outputs
acc = c()
prec = c()
rec = c()
F1 = c()
AUC = c()
predictions = data.frame('ID' = rownames(dataset), 'SFARI' = dataset$SFARI, 'prob' = 0, 'pred' = 0, 'n' = 0)
coefs = data.frame('var' = c('Intercept', colnames(dataset[,-ncol(dataset)])), 'coef' = 0)

for(seed in seeds){
  
  # Run model
  model_output = run_model(p, seed)
  
  # Update outputs
  acc = c(acc, model_output[['acc']])
  prec = c(prec, model_output[['prec']])
  rec = c(rec, model_output[['rec']])
  F1 = c(F1, model_output[['F1']])
  AUC = c(AUC, model_output[['AUC']])
  preds = model_output[['preds']]
  coefs$coef = coefs$coef + model_output[['coefs']]
  update_preds = preds %>% dplyr::select(-ID) %>% mutate(n=1)
  predictions[predictions$ID %in% preds$ID, c('prob','pred','n')] = predictions[predictions$ID %in% 
                                                                      preds$ID, c('prob','pred','n')] +
                                                                    update_preds
}

coefs = coefs %>% mutate(coef = coef/n_iter)
predictions = predictions %>% mutate(prob = prob/n, pred_count = pred, pred = prob>0.5)


rm(p, seeds, update_preds, create_train_test_sets, run_model)


To summarise in a single value the predictions of the models:

test_set = predictions %>% filter(n>0) %>% 
           left_join(dataset %>% mutate(ID = rownames(.)) %>% dplyr::select(ID, GS, MTcor), by = 'ID')
rownames(test_set) = predictions$ID[predictions$n>0]


Performance metrics


Confusion matrix

conf_mat = test_set %>% apply_labels(SFARI = 'Actual Labels', 
                                     prob = 'Assigned Probability', 
                                     pred = 'Label Prediction')

cro(conf_mat$SFARI, list(conf_mat$pred, total()))
 Label Prediction     #Total 
 FALSE   TRUE   
 Actual Labels 
   FALSE  5610 1584   7194
   TRUE  224 209   433
   #Total cases  5834 1793   7627
rm(conf_mat)


Accuracy: Mean = 0.7605 SD = 0.012


Precision: Mean = 0.1131 SD = 0.0099


Recall: Mean = 0.4819 SD = 0.0401


F1 score: Mean = 0.1831 SD = 0.0154


ROC Curve: Mean = 0.6849 SD = 0.0222

pred_ROCR = prediction(test_set$prob, test_set$SFARI)

roc_ROCR = performance(pred_ROCR, measure='tpr', x.measure='fpr')
auc = performance(pred_ROCR, measure='auc')@y.values[[1]]

plot(roc_ROCR, main=paste0('ROC curve (AUC=',round(auc,2),')'), col='#009999')
abline(a=0, b=1, col='#666666')


Lift Curve

lift_ROCR = performance(pred_ROCR, measure='lift', x.measure='rpp')
plot(lift_ROCR, main='Lift curve', col='#86b300')

rm(pred_ROCR, roc_ROCR, AUC, lift_ROCR)




Coefficients


MTcor has a very small coefficient, Gene Significance has a negative coefficient and absGS a positive one

gene_corr_info = dataset %>% mutate('ID' = rownames(dataset)) %>% dplyr::select(ID, MTcor, SFARI) %>% 
                 left_join(assigned_module, by ='ID') %>% mutate(Module = gsub('#','',Module))

coef_info = coefs %>% mutate('feature' = gsub('MM.','',var)) %>% 
            left_join(gene_corr_info, by = c('feature' = 'Module')) %>% 
            dplyr::select(feature, coef, MTcor, SFARI) %>% group_by(feature, coef, MTcor) %>% 
            summarise('SFARI_perc' = mean(SFARI)) %>% arrange(desc(coef))

coef_info %>% dplyr::select(feature, coef) %>% filter(feature %in% c('Intercept','GS','absGS','MTcor')) %>%
              dplyr::rename('Feature' = feature, 'Coefficient' = coef) %>% 
              kable(align = 'cc', 
                    caption = 'Regression Coefficients (filtering MM features)') %>% 
              kable_styling(full_width = F)
Regression Coefficients (filtering MM features)
Feature Coefficient
absGS 0.4405582
MTcor -0.0229410
GS -0.2863086
Intercept -0.6372805


There is a positive relation between the coefficient assigned to the membership of each module and the enrichment (using ORA) in SFARI genes that are assigned to that module

load('./../Data/ORA.RData')

enrichment_SFARI_info = data.frame('Module'=as.character(), 'SFARI_enrichment'=as.numeric())
for(m in names(enrichment_SFARI)){
  m_info = enrichment_SFARI[[m]]
  enrichment = 1-ifelse('SFARI' %in% m_info$ID, m_info$pvalue[m_info$ID=='SFARI'],1)
  enrichment_SFARI_info = enrichment_SFARI_info %>% 
                          add_row(Module = gsub('#','',m), SFARI_enrichment = enrichment)
}

plot_data = coef_info %>% dplyr::rename('Module' = feature) %>% 
            left_join(enrichment_SFARI_info, by = 'Module') %>% filter(!is.na(MTcor))

ggplotly(plot_data %>% ggplot(aes(coef, SFARI_enrichment)) + 
         geom_smooth(method = 'lm', color = 'gray', alpha = 0.1) + 
         geom_point(aes(id = Module), color = paste0('#',plot_data$Module), alpha=0.7) + 
         theme_minimal() + xlab('Coefficient') + 
         ylab('SFARI Genes Enrichment'))
rm(enrichment_old_SFARI, enrichment_DGN, enrichment_DO, enrichment_GO, enrichment_KEGG, enrichment_Reactome,
   m, m_info, enrichment)


There doesn’t seem to be a relation between the coefficient and the correlation of the module and the diagnosis.

This is not a surprise since we knew that there was a negative relation between SFARI genes and Module-Diagnosis correlation from Preprocessing/Gandal/AllRegions/RMarkdowns/20_04_03_WGCNA_modules_EA.html. The fact that there is no relation between coefficient and Module-Diagnosis correlation could even be a good sign that the model is picking some biological signal as well as the SFARI patterns (since the relation with the biological signals is positive)

ggplotly(coef_info %>% dplyr::rename('Module' = feature) %>% filter(!is.na(MTcor)) %>%
              ggplot(aes(coef, MTcor)) +  geom_smooth(method = 'lm', color = 'gray', alpha = 0.1) + 
              geom_point(aes(id = Module), color = paste0('#',coef_info$feature[!is.na(coef_info$MTcor)]), 
                         alpha = 0.7) + 
              theme_minimal() + xlab('Coefficient') + 
              ylab('Module-Diagnosis correlation'))




Analyse Results


Score distribution by SFARI Label


SFARI genes have a higher score distribution than the rest, but the overlap is large

plot_data = test_set %>% dplyr::select(prob, SFARI)

ggplotly(plot_data %>% ggplot(aes(prob, fill=SFARI, color=SFARI)) + geom_density(alpha=0.3) + xlab('Score') +
         geom_vline(xintercept = mean(plot_data$prob[plot_data$SFARI]), color = '#00C0C2', linetype = 'dashed') +
         geom_vline(xintercept = mean(plot_data$prob[!plot_data$SFARI]), color = '#FF7371', linetype = 'dashed') +
         theme_minimal() + ggtitle('Model score distribution by SFARI Label'))


Score distribution by SFARI Score


Even though we didn’t use the actual SFARI Scores to train the model, but instead we grouped them all together, there seems to be a statistically significant positive relation between the SFARI scores and the probability assigned by the model

plot_data = test_set %>% mutate(ID=rownames(test_set)) %>% dplyr::select(ID, prob) %>%
            left_join(original_dataset %>% mutate(ID=rownames(original_dataset)), by='ID') %>%
            mutate(gene.score = ifelse(gene.score=='None', ifelse(ID %in% GO_neuronal$ID, 'Neuronal', 'Others'), 
                                       gene.score)) %>%
            dplyr::select(ID, prob, gene.score) %>% apply_labels(gene.score='SFARI Gene score')

cro(plot_data$gene.score)
 #Total 
 SFARI Gene score 
   1  95
   2  115
   3  223
   Neuronal  539
   Others  6655
   #Total cases  7627
mean_vals = plot_data %>% group_by(gene.score) %>% summarise(mean_prob = mean(prob))

comparisons = list(c('1','2'), c('2','3'), c('3','Neuronal'), c('Neuronal','Others'),
                   c('1','3'), c('3','Others'), c('2','Neuronal'),
                   c('1','Neuronal'), c('2','Others'), c('1','Others'))
increase = 0.08
base = 0.9
pos_y_comparisons = c(rep(base, 4), rep(base + increase, 2), base + 2:5*increase)

plot_data %>% ggplot(aes(gene.score, prob, fill=gene.score)) + 
              geom_boxplot(outlier.colour='#cccccc', outlier.shape='o', outlier.size=3) +
              stat_compare_means(comparisons = comparisons, label = 'p.signif', method = 't.test', 
                                 method.args = list(var.equal = FALSE), label.y = pos_y_comparisons, 
                                 tip.length = .02) +
              scale_fill_manual(values=SFARI_colour_hue(r=c(1:3,8,7))) + 
              ggtitle('Distribution of probabilities by SFARI score') +
              xlab('SFARI score') + ylab('Probability') + theme_minimal() + theme(legend.position = 'none')

rm(mean_vals, increase, base, pos_y_comparisons)


Genes with the highest Probabilities


  • Considering the class imbalance in the test set (1:19), there are many more SFARI scores in here (1:5)

  • High concentration of genes with a SFARI Score of 1

test_set %>% dplyr::select(prob, SFARI) %>% mutate(ID = rownames(test_set)) %>% 
             arrange(desc(prob)) %>% top_n(50, wt=prob) %>%
             left_join(original_dataset %>% mutate(ID=rownames(original_dataset)), by='ID')  %>%
             mutate(gene.score = ifelse(gene.score=='None', ifelse(ID %in% GO_neuronal$ID, 'Neuronal', 'Others'), 
                                        gene.score)) %>%
             left_join(gene_names, by = c('ID'='ensembl_gene_id')) %>%
             dplyr::rename('GeneSymbol' = external_gene_id, 'Probability' = prob, 'ModuleDiagnosis_corr' =MTcor,
                           'GeneSignificance' = GS) %>%
             mutate(ModuleDiagnosis_corr = round(ModuleDiagnosis_corr, 4), Probability = round(Probability, 4), 
                    GeneSignificance = round(GeneSignificance, 4)) %>%
             dplyr::select(GeneSymbol, GeneSignificance, ModuleDiagnosis_corr, Module, Probability,
                           gene.score) %>%
             kable(caption = 'Genes with highest model probabilities from the test set') %>% 
             kable_styling(full_width = F)
Genes with highest model probabilities from the test set
GeneSymbol GeneSignificance ModuleDiagnosis_corr Module Probability gene.score
NCOR2 -0.3804 -0.5674 #F57962 0.8571 Others
BSN -0.1125 -0.1508 #40B500 0.8276 Neuronal
PHACTR4 0.4146 0.7500 #6EB000 0.8226 Others
KMT2D -0.0219 -0.1905 #00C19C 0.8166 Others
ASH1L 0.2512 0.3460 #00B933 0.8132 1
GPATCH2L 0.3333 0.3460 #00B933 0.8131 Others
ARHGAP32 0.2334 0.1879 #FD61D0 0.8125 3
SYNGAP1 -0.2678 -0.6976 #FF64B4 0.8061 1
PLXNA4 -0.0263 0.1879 #FD61D0 0.8053 2
CDK5R1 -0.4394 -0.6976 #FF64B4 0.8038 Neuronal
MYCBP2 0.2898 0.1879 #FD61D0 0.8036 Neuronal
SAMD12 0.1065 0.1879 #FD61D0 0.8020 Others
AMBRA1 -0.0390 -0.1508 #40B500 0.8003 Neuronal
MEGF8 -0.3370 0.0026 #E58700 0.7992 Others
SRRM4 0.1119 0.1879 #FD61D0 0.7947 Others
TROVE2 0.6397 0.7500 #6EB000 0.7930 Others
KMT2C 0.1929 0.2419 #C49A00 0.7929 1
NRXN2 -0.1113 0.0026 #E58700 0.7922 1
SESTD1 -0.0110 0.3460 #00B933 0.7899 Others
PPP1R9A 0.1708 0.3460 #00B933 0.7837 Neuronal
MAP4K4 0.0292 0.3460 #00B933 0.7790 Others
NRXN1 0.0743 0.1879 #FD61D0 0.7763 1
LMTK2 -0.4932 -0.6976 #FF64B4 0.7756 Others
KIF21A 0.1280 0.3460 #00B933 0.7749 Others
HOMER1 0.1537 0.7500 #6EB000 0.7744 3
PHF12 0.2488 0.7335 #FF67A4 0.7737 Others
GSK3B 0.0139 -0.1905 #00C19C 0.7734 Neuronal
PCLO 0.2425 0.3460 #00B933 0.7710 Others
ZFYVE28 -0.2804 -0.5103 #00C083 0.7703 Others
RNF165 0.0710 0.0026 #E58700 0.7701 Others
KMT2A 0.4900 0.2811 #00BD66 0.7695 1
TCF4 0.5478 0.3460 #00B933 0.7693 1
NRXN3 0.2184 0.1879 #FD61D0 0.7685 1
DYNC1H1 -0.1332 -0.1508 #40B500 0.7682 2
BAZ2B 0.3256 0.3460 #00B933 0.7682 1
SLC8A2 -0.4859 -0.6976 #FF64B4 0.7658 Others
FASN -0.4792 -0.6976 #FF64B4 0.7648 Others
TP53INP2 -0.1228 0.1192 #8B93FF 0.7644 Others
PRKCB 0.2515 0.6651 #E86BF2 0.7644 2
REXO1 -0.4670 -0.6022 #998EFF 0.7642 Others
ZDBF2 -0.2698 0.1879 #FD61D0 0.7640 Others
TMEM245 0.4730 0.7500 #6EB000 0.7615 Others
MEF2D -0.3326 -0.1167 #00BC50 0.7613 Others
KIAA1244 0.1185 0.1192 #8B93FF 0.7606 Others
KIF3C -0.2740 -0.5004 #53B400 0.7604 Others
SRGAP3 0.1436 0.7500 #6EB000 0.7599 3
WBSCR17 -0.2650 -0.4648 #7B97FF 0.7598 Others
MINK1 -0.3790 -0.1905 #00C19C 0.7595 Others
NCAM1 0.1477 -0.1508 #40B500 0.7589 Neuronal
PCDH9 0.2477 0.3460 #00B933 0.7588 3





Negative samples distribution


The objective of this model is to identify candidate SFARI genes. For this, we are going to focus on the negative samples (the non-SFARI genes)

negative_set = test_set %>% filter(!SFARI)

negative_set_table = negative_set %>% apply_labels(prob = 'Assigned Probability', 
                                                   pred = 'Label Prediction')

cro(negative_set_table$pred)
 #Total 
 Label Prediction 
   FALSE  5610
   TRUE  1584
   #Total cases  7194

1584 genes are predicted as ASD-related


negative_set %>% ggplot(aes(prob)) + geom_density(color='#F8766D', fill='#F8766D', alpha=0.5) +
                 geom_vline(xintercept=0.5, color='#333333', linetype='dotted') + xlab('Probability') +
                 ggtitle('Probability distribution of the Negative samples in the Test Set') + 
                 theme_minimal()




Probability and Gene Significance


  • There’s a lot of noise, but the probability the model assigns to each gene seems to have a positive relation with the Gene Significance (under-expressed genes having on average the lower probabilities and over-expressed genes the highest) (this pattern was the opposite in Gandal’s dataset)
negative_set %>% ggplot(aes(prob, GS, color = MTcor)) + geom_point() + 
                 geom_smooth(method = 'loess', color = '#666666') +
                 geom_hline(yintercept = 0, color='gray', linetype='dashed') + 
                 xlab('Probability') + ylab('Gene Significance') +
                 scale_color_gradientn(colours=c('#F8766D','white','#00BFC4')) + 
                 ggtitle('Relation between Probability and Gene Significance') + theme_minimal()




Probability and Module-Diagnosis correlation


  • There’s not a strong relation between the Module-Diagnosis correlation of the genes assigned module and the probability assigned by the model

  • The model seems to assign slightly higher probabilities to genes belonging the modules with high module-Dianosis correlations (both positive and negative ones) than to genes belonging to modules with low correlations

negative_set %>% ggplot(aes(MTcor, prob, color=GS)) + geom_point() + 
                 geom_smooth(method='loess', color='#666666') + 
                 geom_hline(yintercept=mean(negative_set$prob), color='gray', linetype='dashed') +
                 scale_color_gradientn(colours=c('#F8766D','#F8766D','white','#00BFC4','#00BFC4')) + 
                 xlab('Modules ordered by their correlation to ASD') + ylab('Model probability') +
                 theme_minimal()




Summarised version, plotting by module instead of by gene


The difference in the trend lines between this plot and the one above is that the one above takes all the points into consideration while this considers each module as an observation by itself, so the top one is strongly affected by big modules and the bottom one treats all modules the same

Again, the model seems to give lower probabilities to genes belonging to modules with a low (absolute) correlation to Diagnosis than to the rest

plot_data = negative_set %>% group_by(MTcor) %>% summarise(mean = mean(prob), sd = sd(prob), n = n()) %>%
            mutate(MTcor_sign = ifelse(MTcor>0, 'Positive', 'Negative')) %>% 
            left_join(original_dataset, by='MTcor') %>%
            dplyr::select(Module, MTcor, MTcor_sign, mean, sd, n) %>% distinct()
colnames(plot_data)[1] = 'ID'

ggplotly(plot_data %>% ggplot(aes(MTcor, mean, size=n, color=MTcor_sign)) + geom_point(aes(id=ID), alpha=0.7) + 
         geom_smooth(method='loess', color='gray', se=FALSE) + geom_smooth(method='lm', se=FALSE) + 
         xlab('Module-Diagnosis correlation') + ylab('Mean Probability by Module') + theme_minimal())




Probability and level of expression


# Gupta dataset
load('./../Data/preprocessed_data.RData')
datExpr = datExpr %>% data.frame
DE_info = DE_info %>% data.frame

There doesn’t seem to be any relation between the level of expression of the genes and the probability assigned to the model despite this happening in Gandal’s dataset

mean_and_sd = data.frame(ID=rownames(datExpr), meanExpr=rowMeans(datExpr), sdExpr=apply(datExpr,1,sd))

plot_data = negative_set %>% left_join(mean_and_sd, by='ID') %>% 
            left_join(original_dataset %>% mutate(ID=rownames(original_dataset)) %>% 
                      dplyr::select(ID, Module), by='ID')
colnames(plot_data)[ncol(plot_data)] = 'Module'

plot_data %>% ggplot(aes(meanExpr, prob)) + geom_point(alpha=0.2, color='#0099cc') + 
              geom_smooth(method='loess', color='gray', alpha=0.3) +  
              geom_smooth(method='lm', color='#808080', se=FALSE, linetype = 'dashed') + 
              xlab('Mean Expression') + ylab('Probability') + 
              ggtitle('Mean expression vs model probability by gene') +
              theme_minimal()

rm(mean_and_sd)

Grouping the genes by module we see there may be a small positive non-linear correlation between mean level of expression and model probability, but it’s still not as strong as in Gandal’s dataset

plot_data2 = plot_data %>% group_by(Module) %>% summarise(meanExpr = mean(meanExpr), meanProb = mean(prob), 
                                                          n=n())

ggplotly(plot_data2 %>% ggplot(aes(meanExpr, meanProb, size=n)) + 
         geom_point(color=plot_data2$Module, alpha = 0.6) + 
         geom_smooth(method='loess', se=TRUE, color='gray', alpha=0.1, size=0.7) + 
         geom_smooth(method='lm', color='#808080', se=FALSE, linetype = 'dashed') + 
         xlab('Mean Level of Expression') + ylab('Average Model Probability') +
         theme_minimal() + theme(legend.position='none') + 
         ggtitle('Mean expression vs model probability by Module'))
rm(plot_data2)




Probability and LFC


There is a relation between probability and LFC, so it IS capturing a bit of true information (because LFC and mean expression were negatively correlated and it still has a positive relation in the model)

  • The relation is stronger in over-expressed in ASD (opposite to the behaviour found in Gandal’s dataset)
plot_data = negative_set %>% left_join(DE_info %>% mutate(ID=rownames(DE_info)), by='ID')

plot_data %>% ggplot(aes(log2FoldChange, prob)) + geom_point(alpha=0.1, color='#0099cc') + 
              geom_smooth(method='loess', color='gray', alpha=0.3) + xlab('LFC') +
              xlab('LFC') + ylab('Probability') +
              theme_minimal() + ggtitle('LFC vs model probability by gene')

  • The relation is stronger in Differentially Expressed genes
p1 = plot_data %>% filter(log2FoldChange<0) %>% mutate(DE = padj<0.05) %>%
     ggplot(aes(log2FoldChange, prob, color=DE)) + geom_point(alpha=0.2) + 
     geom_smooth(method='loess', alpha=0.1) + xlab('') + ylab('Probability') + 
     ylim(c(min(plot_data$prob), max(plot_data$prob))) + 
     theme_minimal() + theme(legend.position = 'none', plot.margin=unit(c(1,-0.3,1,1), 'cm'))

p2 = plot_data %>% filter(log2FoldChange>=0) %>% mutate(DE = padj<0.05) %>% 
     ggplot(aes(log2FoldChange, prob, color=DE)) + geom_point(alpha=0.2) + 
     geom_smooth(method='loess', alpha=0.1) + xlab('') + ylab('Probability') + ylab('') +
     scale_y_continuous(position = 'right', limits = c(min(plot_data$prob), max(plot_data$prob))) +
     theme_minimal() + theme(plot.margin = unit(c(1,1,1,-0.3), 'cm'), axis.ticks.y = element_blank())

grid.arrange(p1, p2, nrow=1, top = 'LFC vs model probability by gene', bottom = 'LFC')

rm(p1, p2)



Conclusion


  • The model doesn’t seem to be capturing the level of expression of the genes as strongly as it did in Gandal’s dataset. This would mean that no bias correction is not as necessary for this dataset, and that the bias in probabilities we observe by SFARI Score is defined by something else.

  • It seems to be capturing true biological signals (based on the GS and the log fold change plots).


Saving results

predictions = test_set %>% left_join(gene_names, by = c('ID' = 'ensembl_gene_id'))

save(predictions, dataset, fit, file='./../Data/Ridge_model.RData')




Session info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] DMwR_0.4.1         doParallel_1.0.15  iterators_1.0.12   foreach_1.5.0     
##  [5] kableExtra_1.1.0   expss_0.10.2       corrplot_0.84      MLmetrics_1.1.1   
##  [9] car_3.0-7          carData_3.0-3      ROCR_1.0-7         gplots_3.0.3      
## [13] caret_6.0-86       lattice_0.20-41    polycor_0.7-10     biomaRt_2.40.5    
## [17] ggpubr_0.2.5       magrittr_1.5       RColorBrewer_1.1-2 gridExtra_2.3     
## [21] viridis_0.5.1      viridisLite_0.3.0  plotly_4.9.2       knitr_1.28        
## [25] forcats_0.5.0      stringr_1.4.0      dplyr_1.0.0        purrr_0.3.4       
## [29] readr_1.3.1        tidyr_1.1.0        tibble_3.0.1       ggplot2_3.3.2     
## [33] tidyverse_1.3.0   
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_1.1.0            RSQLite_2.2.0              
##   [3] AnnotationDbi_1.46.1        htmlwidgets_1.5.1          
##   [5] BiocParallel_1.18.1         pROC_1.16.2                
##   [7] munsell_0.5.0               codetools_0.2-16           
##   [9] miniUI_0.1.1.1              withr_2.2.0                
##  [11] colorspace_1.4-1            Biobase_2.44.0             
##  [13] highr_0.8                   rstudioapi_0.11            
##  [15] stats4_3.6.3                ggsignif_0.6.0             
##  [17] TTR_0.23-6                  labeling_0.3               
##  [19] GenomeInfoDbData_1.2.1      bit64_0.9-7                
##  [21] farver_2.0.3                vctrs_0.3.1                
##  [23] generics_0.0.2              ipred_0.9-9                
##  [25] xfun_0.12                   R6_2.4.1                   
##  [27] GenomeInfoDb_1.20.0         locfit_1.5-9.4             
##  [29] bitops_1.0-6                DelayedArray_0.10.0        
##  [31] assertthat_0.2.1            promises_1.1.0             
##  [33] scales_1.1.1                nnet_7.3-14                
##  [35] ggExtra_0.9                 gtable_0.3.0               
##  [37] timeDate_3043.102           rlang_0.4.6                
##  [39] genefilter_1.66.0           splines_3.6.3              
##  [41] lazyeval_0.2.2              ModelMetrics_1.2.2.2       
##  [43] acepack_1.4.1               broom_0.5.5                
##  [45] checkmate_2.0.0             yaml_2.2.1                 
##  [47] reshape2_1.4.4              abind_1.4-5                
##  [49] modelr_0.1.6                crosstalk_1.1.0.1          
##  [51] backports_1.1.8             quantmod_0.4.17            
##  [53] httpuv_1.5.2                Hmisc_4.4-0                
##  [55] tools_3.6.3                 lava_1.6.7                 
##  [57] ellipsis_0.3.1              BiocGenerics_0.30.0        
##  [59] Rcpp_1.0.4.6                plyr_1.8.6                 
##  [61] base64enc_0.1-3             progress_1.2.2             
##  [63] zlibbioc_1.30.0             RCurl_1.98-1.2             
##  [65] prettyunits_1.1.1           rpart_4.1-15               
##  [67] zoo_1.8-8                   S4Vectors_0.22.1           
##  [69] SummarizedExperiment_1.14.1 haven_2.2.0                
##  [71] cluster_2.1.0               fs_1.4.0                   
##  [73] data.table_1.12.8           openxlsx_4.1.4             
##  [75] reprex_0.3.0                matrixStats_0.56.0         
##  [77] hms_0.5.3                   mime_0.9                   
##  [79] evaluate_0.14               xtable_1.8-4               
##  [81] XML_3.99-0.3                rio_0.5.16                 
##  [83] jpeg_0.1-8.1                readxl_1.3.1               
##  [85] shape_1.4.4                 IRanges_2.18.3             
##  [87] compiler_3.6.3              KernSmooth_2.23-17         
##  [89] crayon_1.3.4                htmltools_0.4.0            
##  [91] mgcv_1.8-31                 later_1.0.0                
##  [93] Formula_1.2-3               geneplotter_1.62.0         
##  [95] lubridate_1.7.4             DBI_1.1.0                  
##  [97] dbplyr_1.4.2                MASS_7.3-51.6              
##  [99] Matrix_1.2-18               cli_2.0.2                  
## [101] gdata_2.18.0                gower_0.2.1                
## [103] GenomicRanges_1.36.1        pkgconfig_2.0.3            
## [105] foreign_0.8-76              recipes_0.1.10             
## [107] xml2_1.2.5                  annotate_1.62.0            
## [109] webshot_0.5.2               XVector_0.24.0             
## [111] prodlim_2019.11.13          rvest_0.3.5                
## [113] digest_0.6.25               rmarkdown_2.1              
## [115] cellranger_1.1.0            htmlTable_1.13.3           
## [117] curl_4.3                    shiny_1.4.0.2              
## [119] gtools_3.8.2                lifecycle_0.2.0            
## [121] nlme_3.1-147                jsonlite_1.7.0             
## [123] fansi_0.4.1                 pillar_1.4.4               
## [125] fastmap_1.0.1               httr_1.4.1                 
## [127] survival_3.1-12             xts_0.12-0                 
## [129] glue_1.4.1                  zip_2.0.4                  
## [131] png_0.1-7                   glmnet_3.0-2               
## [133] bit_1.1-15.2                class_7.3-17               
## [135] stringi_1.4.6               blob_1.2.1                 
## [137] DESeq2_1.24.0               latticeExtra_0.6-29        
## [139] caTools_1.18.0              memoise_1.1.0